#===============================================================================
# This code is for replicating the key results on the analyses of robusta coffee
# response to climate stress under early and late flowering conditions. The key results
# and code are provided here. For further details please contact 

# Dr Jarrod Kath, jarrod.kath@usq.edu.au,
# Centre for Applied Climate Sciences, University of Southern Queensland
#
#===============================================================================

###########################################
##1. Load in the dataset
##########################################

setwd("....")

load("Robusta_flower_df.Rdata")

###Give the dataset a shorter name for convenience 
viet<-Robusta_flower_df


################
##2. Create early and late flowering dataframes
#################

viet_early<-subset(viet, f_doy_anom.x<0)

viet_late<-subset(viet, f_doy_anom.x>0)


#####################################################
##3. Transform the different responses as SEM can be sensitive to normality
#######################################

library(bestNormalize)##



###All / overall model transformation
(orderNorm_obj_fdoy_Viet <- orderNorm(viet$f_doy_anom.x))
(orderNorm_obj_gs_Viet <- orderNorm(viet$gs_anom.x))
(orderNorm_obj_yld_Viet <- orderNorm(viet$yld.x))


viet$trans_fdoy<-orderNorm_obj_fdoy_Viet$x.t
viet$trans_gsdoy<-orderNorm_obj_gs_Viet$x.t
viet$trans_yld<-orderNorm_obj_yld_Viet$x.t



###Early model transformation

(orderNorm_obj_fdoy_Viet_early <- orderNorm(viet_early$f_doy_anom.x))
(orderNorm_obj_gs_Viet_early <- orderNorm(viet_early$gs_anom.x))
(orderNorm_obj_yld_Viet_early <- orderNorm(viet_early$yld.x))

##

viet_early$trans_fdoy<-orderNorm_obj_fdoy_Viet_early$x.t
viet_early$trans_gsdoy<-orderNorm_obj_gs_Viet_early$x.t
viet_early$trans_yld<-orderNorm_obj_yld_Viet_early$x.t
# 
# 

###Late model transformation
(orderNorm_obj_fdoy_Viet_late <- orderNorm(viet_late$f_doy_anom.x))
(orderNorm_obj_gs_Viet_late <- orderNorm(viet_late$gs_anom.x))
(orderNorm_obj_yld_Viet_late <- orderNorm(viet_late$yld.x))

##

viet_late$trans_fdoy<-orderNorm_obj_fdoy_Viet_late$x.t
viet_late$trans_gsdoy<-orderNorm_obj_gs_Viet_late$x.t
viet_late$trans_yld<-orderNorm_obj_yld_Viet_late$x.t
# 






#######################################
##3. Centre variables so can easily compare effect size after Gelmand and Hill 2006
#######################################

###NOTE throughout climate variables with ".x" at the end designate those calculated from the flowering period,
### those with ".y" at the end those from the fruiting period and those without and ".x" or ".y" are climate variables from 
### the growing season. Please see Table 1 in the paper for the definition of these phases.


##Centre variables
cent<- function(x) {((x) - mean(x, na.rm=T))/(1*sd(x, na.rm=T))}
reversecent<-function(x,y) {((x)*1*sd(y, na.rm=T)) + mean(y, na.rm=T)}



##All
viet$cent_tmin30_rollmean <-cent(viet$tmin30_rollmean)
viet$cent_tmax30_rollmean <-cent(viet$tmax30_rollmean)
viet$cent_rain30_rollsum <-cent(viet$rain30_rollsum)


viet$cent_tmin30_rollmean.x <-cent(viet$tmin30_rollmean.x)
viet$cent_tmax30_rollmean.x <-cent(viet$tmax30_rollmean.x)
viet$cent_rain30_rollsum.x <-cent(viet$rain30_rollsum.x)

viet$cent_tmin30_rollmean.y <-cent(viet$tmin30_rollmean.y)
viet$cent_tmax30_rollmean.y <-cent(viet$tmax30_rollmean.y)
viet$cent_rain30_rollsum.y <-cent(viet$rain30_rollsum.y)

viet$cent_tmings_rollmean <-cent(viet$tmings_rollmean)
viet$cent_tmaxgs_rollmean <-cent(viet$tmaxgs_rollmean)
viet$cent_raings_rollsum <-cent(viet$raings_rollsum)

viet$cent_irrig<-cent(viet$Tot_irrig.x)
viet$cent_fert<-cent(viet$Tot_fertilizer)
viet$cent_age<-cent(viet$Age_plant)



###early
viet_early$cent_tmin30_rollmean <-cent(viet_early$tmin30_rollmean)
viet_early$cent_tmax30_rollmean <-cent(viet_early$tmax30_rollmean)
viet_early$cent_rain30_rollsum <-cent(viet_early$rain30_rollsum)


viet_early$cent_tmin30_rollmean.x <-cent(viet_early$tmin30_rollmean.x)
viet_early$cent_tmax30_rollmean.x <-cent(viet_early$tmax30_rollmean.x)
viet_early$cent_rain30_rollsum.x <-cent(viet_early$rain30_rollsum.x)

viet_early$cent_tmin30_rollmean.y <-cent(viet_early$tmin30_rollmean.y)
viet_early$cent_tmax30_rollmean.y <-cent(viet_early$tmax30_rollmean.y)
viet_early$cent_rain30_rollsum.y <-cent(viet_early$rain30_rollsum.y)

viet_early$cent_tmings_rollmean <-cent(viet_early$tmings_rollmean)
viet_early$cent_tmaxgs_rollmean <-cent(viet_early$tmaxgs_rollmean)
viet_early$cent_raings_rollsum <-cent(viet_early$raings_rollsum)

viet_early$cent_irrig<-cent(viet_early$Tot_irrig.x)
viet_early$cent_fert<-cent(viet_early$Tot_fertilizer)
viet_early$cent_age<-cent(viet_early$Age_plant)


##late
viet_late$cent_tmin30_rollmean <-cent(viet_late$tmin30_rollmean)
viet_late$cent_tmax30_rollmean <-cent(viet_late$tmax30_rollmean)
viet_late$cent_rain30_rollsum <-cent(viet_late$rain30_rollsum)


viet_late$cent_tmin30_rollmean.x <-cent(viet_late$tmin30_rollmean.x)
viet_late$cent_tmax30_rollmean.x <-cent(viet_late$tmax30_rollmean.x)
viet_late$cent_rain30_rollsum.x <-cent(viet_late$rain30_rollsum.x)

viet_late$cent_tmin30_rollmean.y <-cent(viet_late$tmin30_rollmean.y)
viet_late$cent_tmax30_rollmean.y <-cent(viet_late$tmax30_rollmean.y)
viet_late$cent_rain30_rollsum.y <-cent(viet_late$rain30_rollsum.y)

viet_late$cent_tmings_rollmean <-cent(viet_late$tmings_rollmean)
viet_late$cent_tmaxgs_rollmean <-cent(viet_late$tmaxgs_rollmean)
viet_late$cent_raings_rollsum <-cent(viet_late$raings_rollsum)

viet_late$cent_irrig<-cent(viet_late$Tot_irrig.x)
viet_late$cent_fert<-cent(viet_late$Tot_fertilizer)
viet_late$cent_age<-cent(viet_late$Age_plant)



##4.  SEM models
#see here https://cougrstats.wordpress.com/2019/12/03/sem-lesson-for-r-working-group/ adding correlated errors
#and here https://cran.r-project.org/web/packages/piecewiseSEM/vignettes/piecewiseSEM.html 

##Just provided overall model here, but same structure used for early, late model


library(piecewiseSEM)
library(nlme)

coffee.list_M1_early <- list(
  
  ####Flowering day anom.  model
  lme(trans_fdoy ~   cent_rain30_rollsum + cent_tmin30_rollmean + cent_tmax30_rollmean + cent_age  + cent_fert + cent_irrig, random = ~ 1 | prov.x /SITE.x, na.action = na.omit,correlation = corCAR1(form = ~year),
      data = viet),
  
  
  ####Growing season anom. model
  lme(trans_gsdoy ~ trans_fdoy + cent_age  + cent_fert + cent_irrig, random = ~ 1 | prov.x /SITE.x, na.action = na.omit,correlation = corCAR1(form = ~year),
      data = viet),

  
  ##Correlated errors
  
  cent_raings_rollsum %~~% trans_fdoy,
  cent_rain30_rollsum.x %~~% trans_fdoy,
  cent_tmaxgs_rollmean %~~% trans_fdoy,
  cent_tmings_rollmean %~~% trans_fdoy,
  cent_tmax30_rollmean.x %~~% trans_fdoy,
  cent_tmin30_rollmean.x %~~% trans_fdoy,

  cent_rain30_rollsum.x %~~% trans_gsdoy,
  cent_tmings_rollmean %~~% trans_gsdoy,
  cent_tmin30_rollmean.x %~~% trans_gsdoy,
  cent_tmax30_rollmean.x %~~% trans_gsdoy,
  cent_tmax30_rollmean %~~% trans_gsdoy,
  cent_tmin30_rollmean %~~% trans_gsdoy,
  cent_rain30_rollsum %~~% trans_gsdoy,
  trans_gsdoy %~~% cent_tmaxgs_rollmean,
  trans_gsdoy %~~% cent_tmings_rollmean,
  trans_gsdoy %~~% cent_raings_rollsum,


  cent_tmax30_rollmean %~~% cent_tmin30_rollmean,
  cent_tmin30_rollmean.x %~~% cent_tmax30_rollmean.x,
  cent_tmaxgs_rollmean %~~% cent_tmings_rollmean,
  cent_tmax30_rollmean %~~% cent_tmax30_rollmean.x,
  cent_tmin30_rollmean %~~% cent_tmin30_rollmean.x,
  cent_rain30_rollsum %~~% cent_rain30_rollsum.x,
  cent_raings_rollsum %~~% cent_rain30_rollsum,
  cent_irrig %~~% cent_rain30_rollsum.x,
  cent_tmax30_rollmean %~~% cent_tmaxgs_rollmean,
  cent_tmin30_rollmean %~~% cent_tmings_rollmean,


  cent_tmax30_rollmean %~~% trans_yld,
  cent_tmin30_rollmean %~~% trans_yld,
  cent_rain30_rollsum %~~% trans_yld,
  

  ####Yield response model
  
  lme(trans_yld ~ trans_fdoy + trans_gsdoy+
        
        cent_tmin30_rollmean.x +
        cent_tmax30_rollmean.x +
        cent_rain30_rollsum.x +
        
        
        cent_rain30_rollsum.x:cent_tmax30_rollmean.x +
        cent_rain30_rollsum.x:cent_tmin30_rollmean.x +

        
        cent_tmaxgs_rollmean +
        cent_tmings_rollmean +
        cent_raings_rollsum +
        
        
        cent_raings_rollsum:cent_tmaxgs_rollmean +
        cent_raings_rollsum:cent_tmings_rollmean +
        
        cent_age  + cent_fert + cent_irrig
  

      ,
      
      
      random = ~ 1 | prov.x /SITE.x, na.action = na.omit,correlation = corCAR1(form = ~year),
      data = viet)
  
  
)



coffee.list_M1_early.psem <- as.psem(coffee.list_M1_early)
(new.summary <- summary(coffee.list_M1_early.psem, .progressBar = T))




############################################
####Main result from paper
##Figure 3
########################################

library(DiagrammeR)
library(DiagrammeRsvg)
library(rsvg)
###Overall model
##maybe just make non-sig variables a dashed grey line

overall_graph<-grViz("
digraph a_nice_graph {

# node definitions with substituted label text
node [fontname = Helvetica]
a [label = '@@1'] [shape=sqaure]
b [label = '@@2'] [shape=sqaure]
c [label = '@@3'] [shape=sqaure]
d [label = '@@4']
e [label = '@@5']

f [label = '@@6']


g [label = '@@7'] [shape=sqaure]
h [label = '@@8'] [shape=sqaure]
i [label = '@@9'] [shape=sqaure]

m [label = '@@10'] [shape=sqaure]
n [label = '@@11'] [shape=sqaure]

j [label = '@@12'] [shape=sqaure]
k [label = '@@13'] [shape=sqaure]
l [label = '@@14'] [shape=sqaure]



o [label = '@@15'] [shape=sqaure]
p [label = '@@16'] [shape=sqaure]


q [label = '@@17'] [shape=diamond]
r [label = '@@18'] [shape=diamond]
s [label = '@@19'] [shape=diamond]


# edge definitions with the node IDs

a -> d [color = grey] [penwidth=1] [minlen=2] [fontsize=14] [fontcolor=black] [style=dashed]
b -> d [color = lightblue] [penwidth=2.74] [label= 0.274 , minlen=2] [fontsize=14] [fontcolor=black]
c -> d [color = red] [penwidth=7.31] [label=-0.731, minlen=2] [fontsize=14] [fontcolor=black]
 
d -> e [color = red] [penwidth=8.12] [label=-0.812, minlen=2] [fontsize=14] [fontcolor=black]

e -> f [color = grey] [penwidth=1] [minlen=2] [fontsize=14] [fontcolor=black] [style=dashed]
d -> f [color = red] [penwidth=0.602] [label=-0.060, minlen=2] [fontsize=14] [fontcolor=black] 


g -> f [color = red] [penwidth=2.366] [label=-0.237, minlen=4] [fontsize=14] [fontcolor=black]
h -> f [color = red] [penwidth=2.14] [label=-0.214, minlen=4] [fontsize=14] [fontcolor=black]
i -> f [color = lightblue] [penwidth=4.51] [label=0.451, minlen=4] [fontsize=14] [fontcolor=black]


j -> f [color = red] [penwidth=1.40] [label=-0.140, minlen=2] [fontsize=14] [fontcolor=black]
k -> f [color = red] [penwidth=8.17] [label=-0.817, minlen=2] [fontsize=14] [fontcolor=black]
l -> f [color = lightblue] [penwidth=1.77] [label=0.177, minlen=2] [fontsize=14] [fontcolor=black]

m -> f [color = red] [penwidth=2.37] [label=-0.237, minlen=3] [fontsize=14] [fontcolor=black]
n -> f [color = lightblue] [penwidth=3.01] [label=0.301, minlen=3] [fontsize=14] [fontcolor=black]


o -> f [color = grey] [penwidth=1] [minlen=1] [fontsize=14] [fontcolor=black] [style=dashed]
p -> f [color = grey] [penwidth=1] [minlen=1] [fontsize=14] [fontcolor=black] [style=dashed]
 

q -> f [color = lightblue] [penwidth=8.32] [label=0.832, minlen=1] [fontsize=14] [fontcolor=black] 
r -> f [color = lightblue] [penwidth=3.841] [label=0.384, minlen=1] [fontsize=14] [fontcolor=black] 
s -> f [color = lightblue] [penwidth=4.43] [label=0.443, minlen=1] [fontsize=14] [fontcolor=black] 
 
q -> d [color = red] [penwidth=1.67] [label=-0.167, minlen=1] [fontsize=14] [fontcolor=black] 
r -> d [color = lightblue] [penwidth=1.07] [label=0.107, minlen=1] [fontsize=14] [fontcolor=black] 
s -> d [color = grey] [penwidth=1] [minlen=1] [fontsize=14] [fontcolor=black][style=dashed] 
 
q -> e [color = lightblue] [penwidth=0.888] [label=0.088, minlen=1] [fontsize=14] [fontcolor=black] 
r -> e [color = grey] [penwidth=1] [minlen=1] [fontsize=14] [fontcolor=black] [style=dashed] 
s -> e [color = grey] [penwidth=1] [minlen=1] [fontsize=14] [fontcolor=black] [style=dashed] 
 
 
 
}

[1]: paste0('Pre-flower\\n ','rain')
[2]: paste0('Pre-flower\\n ','tmax')
[3]: paste0('Pre-flower\\n ','tmin')
[4]: paste0('Flowering day\\n ','anomaly')
[5]: paste0('Growing season\\n ','length anomaly')
[6]: paste0('Yield')
[7]: paste0('Flower\\n ','rain')
[8]: paste0('Flower\\n ','tmax')
[9]: paste0('Flower\\n ','tmin')
[10]: paste0('Flower rain X\\n ','Flower tmax')
[11]: paste0('Flower rain X\\n ','Flower tmin')
[12]: paste0('Growing\\n ','rain')
[13]: paste0('Growing\\n ','tmax')
[14]: paste0('Growing\\n ','tmin')
[15]: paste0('Growing rainX\\n ','Growing tmax')
[16]: paste0('Growing rain X\\n ','Growing tmin')
[17]: paste0('Tree\\n ','age')
[18]: paste0('Total\\n ','irrigation')
[19]: paste0('Total\\n ','fertilizer')


")

overall_graph



# 2. Convert to SVG, then save as png
tmp = DiagrammeRsvg::export_svg(overall_graph)
tmp = charToRaw(tmp) # flatten
rsvg::rsvg_png(tmp, "overall_graph_v3.tiff") # saved graph as png in current working directory









